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We present the results of a large scale computer simulation of su- 
percooled silica. We find that at high temperatures the diffusion 
constants show a non-Arrhenius temperature dependence whereas at 
low temperature this dependence is also compatible with an Arrhe- 
nius law. We demonstrate that at low temperatures the intermediate 
scattering function shows a two-step relaxation behavior and that it 
obeys the time temperature superposition principle. We also discuss 
the wave- vector dependence of the nonergodicity parameter and the 
time and temperature dependence of the non-Gaussian parameter. 



1 Introduction 

In the last few years ample evidence was given that computer simulations are 
a very useful method to study the dynamics of supercooled liquids (Barrat and 
Klein 1991, Yip 1995, Kob 1995, Paul and Baschnagel 1995). Apart from some rare 
exceptions, such as, e.g., the recent simulation of supercooled water by Sciortino et 
al. (Gallo et al. 1996, Sciortino et al. 1996), most of these simulations were done 
for liquids for which the local structure is similar to the one of a closed packed 
hard sphere system. Thus the dynamics of systems in which the particles form 
an open network structure, such as Si02 or B2O3, has hardly been investigated 
with computer simulations at all, although there is a large body of experiments in 
which such materials have been studied (see, e.g., the papers in Ngai 1994). In this 
paper we present some of the results of a large scale molecular dynamics computer 
simulation in which we investigated the dynamics of Si02 in its supercooled state. 

2 Model 

The model we use to describe silica is given by the potential proposed a few years 
ago by van Beest et al. (van Beest, Kramer and van Santen 1990), and is of the 
form 

' ij ' ij 

The values of the various parameters can be found in the original publication. This 
potential has been shown to give a correct description of the different crystalline 
phases of silica (Tse and Klug 1991; Tse, Klug and Allan 1995) and in a recent 
publication Vollmayr et al. (Vollmayr, Kob and Binder 1996) have shown that it 
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is also able to reproduce many static properties of amorphous silica and that this 
system forms a network of tetrahedra, each of which has a silicon atom in its center 
and which are interconnected by bridging oxygen atoms. Thus we conclude that 
this model gives a quite realistic description of this network glassformer. Note that 
in the simulation of Vollmayr et al. the non-Coulombic interactions were truncated 
and shifted at a distance of 5.5A and in the present work we have done so likewise. 

The simulations were done at constant volume and the density of the system 
was fixed to 2.3 g/cm 3 . In order to avoid finite size effect in the dynamics (Horbach, 
Kob, Binder and Angell 1996), the size of the system was relatively large and 
consisted of 8016 ions. The equations of motion were integrated with the velocity 
form of the Verlet algorithm and the Coulombic contributions to the potential 
were calculated with the Ewald summation. The time step of the integration was 
1.6 fs and the longest runs were 4- 10 6 time steps, thus 6.4 ns. At each temperature 
the system was equilibrated with a stochastic heat bath for a time which was at 
least as long as the subsequent production run. The temperatures investigated 
were 6100 K, 5200 K, 4700 K, 4300 K, 4000 K, 3900 K, 3760 K, 3580 K, 3400 K, 
3250 K, 3100 K, 3000 K and 2900 K, the lowest temperature at which the system 
could be equilibrated within the time span of our simulation. More details on the 
simulation will be given elsewhere (Horbach and Kob 1997). 

3 Results 

One of the simplest quantities to describe the dynamics of the system is the tracer 
diffusion constant D which can be computed from the mean square displacement 
of a tagged particle. In Fig. [I] we show D for the silicon and oxygen atoms in an 
Arrhenius plot. We recognize that at high and intermediate temperatures the diffu- 
sion constants clearly show a deviation from the Arrhenius behavior which is found 
experimentally, although at lower temperatures than we investigated (Brebec et 
al. 1980, Mikkelsen 1984). However, at the lowest temperature considered here, 
our data is also compatible with an Arrhenius temperature dependence (bold solid 
lines). The activation energies obatined here, 4.45 eV for O and 4.9 eV for Si, 
are in surprisingly good agreement with the experimental values (4.7 eV for O 
(Mikkelsen 1984) and 6 eV for Si (Brebec et al. 1980)). 

Motivated by the proposition of the so-called mode-coupling theory (MCT) 
(Gotze 1991, Gotze and Sjogren 1992), we also attempted to fit our data for inter- 
mediate and low temperatures with a power-law, i.e. D oc (T — T c ) 7 . The result 
of this fit is included in the figure as well (thin solid lines) and we recognize that 
this functional form is able to describe the data very well. Note that the critical 
temperature T c = 2260 K is the same for both types of atoms. (As one might 
have expected, for strong glass formers, T c is well above the experimental glass 
transition temperature T g « 1450 K). Hence we can infer that at our lowest tem- 
peratures the diffusion constants make a crossover from a power-law behavior to 
an Arrhenius behavior. Thus from our simulation we can make the conjecture that 
at high enough temperatures real silica will make a transition from being a strong 
glass former (Angell 1985), which shows an Arrhenius temperature dependence 
of the transport coefficients, to a fragile glass former, which shows a power-law 
like temperature dependence. This prediction is also in agreement with recent 
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propositions of Rossler and Sokolov (Rossler and Sokolov 1996). 

In Fig. |2| we show the time dependence of the incoherent intermediate scatter- 
ing function F s (q,t) for silicon for all temperatures investigated. The wave- vector 
is 1.7 A -1 , the location of the first sharp diffraction peak in the structure factor. 
We see, Fig. 0a, that even at relatively high temperatures F s (q, t) shows a shoulder 
at around 0.6 ps, which develops into a plateau upon further cooling. Thus we find 
that also this network glassformer shows the typical two-step relaxation behavior 
that is known for fragile glasses and whose existence is one of the predictions of 
MCT. For the lowest temperatures we also observe that the correlation function 
shows a dip around 0.2 ps. This feature is likely related to the so-called boson 
peak, which has also been observed in n-scattering experiments of vitreous silica 
(Buchenau et al. 1986). 

We define the a-relaxation time r(T) as the time it takes the correlation 
function to decay to e _1 of its initial value. MCT predicts that if a time correlation 
function is plotted versus the rescaled time t/r(T), one obtains a master curve 
[time-temperature superposition principle (TTSP)]. That this is indeed the case is 
demonstrated in Fig. ^|b, where we show the same curves as in Fig. ^|a, but versus 
t/r(T). For rescaled times t/r > 1, i.e. in the a-relaxation regime, the curves for 
the different temperatures fall nicely on top of each other, thus showing that in 
this time regime the TTSP holds. This is not quite the case for t/r < 1, i.e. in the 
/3-relaxation region. However, a closer inspection of the individual curves shows 
that the reason for the lack of a perfect scaling is likely the fact that the curves 
show the aforementioned dip at early times. At the very lowest temperatures, the 
location of this dip has been moved to such small rescaled times, that the scaling 
of the curves works very well again. Thus we find that at low enough temperatures 
the TTSP works also in the /3-relaxation regime. Similar conclusions hold also for 
other wave- vectors and the oxygen atoms (Horbach and Kob 1997). 

The intermediate scattering function can be measured relatively easily in scat- 
tering experiments. This is not the case for the so-called non-Gaussian parame- 
ters aa(t), which are a measure for how strongly F s (q,t) deviates from a Gaus- 
sian behavior (Rahman 1964), whereas these functions can be calculated easily in 
a computer simulation. (In recent n-scattering experiments by Buchenau et al. 
and Zorn the time dependence of a 2 was determined (Buchenau et al 1996, Zorn 
1996)). In Fig. |]we show the time dependence of a 2 = 3(r 4 (t))/5(r 2 (t)) 2 — 1 and 
a 3 = 3(r 6 (t))/35(r 2 (t)) 3 — 1 for the oxygen atoms, were r(t) is the probability that 
a tagged particle travels a distance r in time t. 

From this figure we see that for short times en is zero, as is should be, and then 
starts to increase when the system enters the /3-relaxation region. After having 
reached a maximum value at a time which corresponds to approximately the end 
of the /5-relaxation regime, the curves start to drop back to zero. The height of this 
maximum increases with decreasing temperature, thus showing that the Gaussian 
approximation breaks down more and more the lower the temperature is. It is 
also interesting that for very short times, m 0.04 ps, a.\i<i show a small peak. A 
similar peak was also found in a computer simulation of a Lennard- Jones system 
(Kob and Andersen 1995a), although with much smaller amplitude. Since this 
latter system is not a network forming glass and quite fragile, we conclude that 
this feature is probably not closely connected to the details of the structure of the 
system but is more likely a feature that is related to the dynamics of supercooled 



liquids. 

The last quantity we investigate is the wave-vector dependence of the noner- 
godicity parameter, i.e. the height of the plateau in a time correlation function at 
low temperatures (Gotze 1991, Gotze and Sjogren 1992). For the case of the coher- 
ent and incoherent intermediate scattering function this is thus nothing else than 
the Debye- Waller and Lamb-Mossbauer factor, respectively. In Fig. |] we show the 
nonergodicity parameters fj^ s \q) and f c (q) for the incoherent scattering function 
of the silicon atoms and the coherent scattering function for the silicon-silicon cor- 
relation, respectively. We see that shows a Gaussian like dependence on q, 
as it has already been observed for the case of simple liquids (Kob and Andersen 
1995b), and is also predicted by MCT. The nonergodicity parameter for the coher- 
ent scattering function shows several maxima and minima, the location of which 
are the same as the ones observed in the structure factor. Also this observation is 
in qualitative agreement with the one made for simple liquids and the prediction 
of MCT. It is also interesting to note that at a wave-vector of approximately 8A _1 
the nonergodicity parameters have decayed to essentially zero, whereas the three 
partial structure factors show a significant wave- vector dependence even for larger 
values of q. This is in contrast to what has been found in the case of simple liq- 
uids (Kob and Andersen 1995b), where the value of the wave- vector at which the 
nonergodicity parameters become zero and the ones at which the partial structure 
factors become constant are almost identical. It would be interesting to under- 
stand whether this observation is just a particularity of silica or whether it is a 
general property of network forming glasses. 
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Figures 



Figure 1: Arrhenius plot of the diffusion constant for the silicon atoms (filled 
squares) and oxygen atoms (open circles). Thin solid lines: Fit with a power-law. 
Bold solid lines: Fit with an Arrhenius-law. 

Figure 2: Incoherent intermediate scattering function for the silicon atoms vs. 
time (a) and rescaled time (b) for all temperatures investigated. 

Figure 3: Time dependence of the non-Gaussian parameters a<i (a) and a 3 (b) for 
all temperatures investigated (main figure). Inset: Same data at short times. 

Figure 4: Wave-vector dependence of the nonergodicity parameter for the inco- 
herent scattering function of the silicon atoms (open circles) and of the coherent 
scattering function of the silicon-silicon correlation (filled squares). 
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